#### https://r2clr.codeplex.com/wikipage?title=Spatial-temporal%20model
rclrRoot = 'f:/codeplex/r2clr' # set this to your own location

library(rClr)

# Use the couple of following lines if you want to run basic unit tests prior to undertaking the tutorial.
# library(testthat)
# test_package('rClr')

library(zoo) # A fantastic package to work with time series

clrLoadAssembly(file.path( rclrRoot, 'doc/src/EnvModellingSample/bin/Release/EnvModellingSample.dll'))
# or
clrLoadAssembly(file.path( rclrRoot, 'doc/src/EnvModellingSample/bin/Debug/EnvModellingSample.dll'))

# Load observed data for a catchment (a.k.a watershed)
load(file=file.path( rclrRoot, 'doc/data/catData.RData'))
rainfall = catData[,1]
pet = catData[,2] # potential evaporation (maybe not quite the right input, but will do for a tutorial)
obsRunoff = catData[,3]
str(rainfall)

simulation = clrCallStatic('EnvModellingSample.SimulationFactory,EnvModellingSample', 'CreateAwbmSimulation')

clrReflect(simulation)

### A few R functions to contain the proliferations of strings.
playInput <- function(simulation, inputName, data) {
  clrCall(simulation, 'Play', inputName, as.numeric(data))
}

recordOutput <- function(simulation, outputName) {
  clrCall(simulation, 'Record', outputName)
}

getRecordedOutput <- function(simulation, outputName) {
  numericVec = clrCall(simulation, 'GetRecorded', outputName)
  zoo(numericVec, index(rainfall))
}

startInd = as.integer(0)
endInd = as.integer(length(rainfall)-1)
clrCall(simulation, 'SetTimeSpan', startInd, endInd)

playInput(simulation, 'Rainfall', rainfall)
playInput(simulation, 'Evapotranspiration', pet)

recordOutput(simulation, 'Runoff')

clrCall(simulation, 'Execute')

model = clrCall(simulation, 'get_TsModel')

## Let's have a look at what the default output of the model looks like, given inputs and current parameters
calcRunoff = getRecordedOutput(simulation, 'Runoff')
plot(merge(rainfall, calcRunoff))

plotWindow = as.Date(c('1990-01-01','1992-12-31'))
calcAndObs = window(merge(obsRunoff, calcRunoff), start=plotWindow[1], end=plotWindow[2])
plot( calcAndObs, plot.type='single', col=c('blue','red'))


### Functions used to summarise the goodness of fit between observations and omodel predictions.

#  http://en.wikipedia.org/wiki/Nash%E2%80%93Sutcliffe_model_efficiency_coefficient
nse <- function(calc, obs)
{
  valid = which(!is.na(obs))
  obsMean = mean(obs,na.rm=TRUE)
  differences = calc-obs
  differences = differences[valid]
  diffAboutMean = obs-obsMean
  diffAboutMean = diffAboutMean[valid]
  1-(sum(differences*differences) / sum(diffAboutMean*diffAboutMean))
}
# tests:
# > nse(1:5, 1:5)
# [1] 1
# > nse(1:5, 1:5*1.1)
# [1] 0.9545455
# > nse(1:5*1.1, c(1,2,NA,4,5))
# [1] 0.954

startObj=as.Date('1985-01-01')
endObj=as.Date('1999-12-31')

objectiveFunc <- function(calc,obs) {
  calcw = window(calc, start=startObj, end=endObj)
  obsw = window(obs, start=startObj, end=endObj)
  nse(calc, obs)
}

getResponse <- function(parameters, simulation, f) {
  for(i in 1:nrow(parameters)) clrCall(simulation, 'SetValue', parameters$name[i], parameters$value[i])
  clrCall(simulation, 'Execute')
  calcRunoff = getRecordedOutput(simulation, 'Runoff')
  f(calcRunoff, obsRunoff)
}


objectiveFunc(calcRunoff, obsRunoff)
# [1] 0.5978418


### Generation of candidate parameter sets

template <- data.frame(name=c('BFI','KSurf','KBase','C1','C2','C3'), min=c(1e-3,1e-3,1e-3,0,0,0), max=c(0.999,0.999,0.999,50,200,500), value=rep(0.5,6), stringsAsFactors=FALSE)

createCandidateParam <- function(template) {
  res <- template
  res$value = runif(nrow(template), template$min, template$max)
  res
}

### Evaluate the goodness of fit for 1000 parameterisations

set.seed(0)
n=1000
# n=10
m = nrow(template)+1
responses <- as.data.frame(matrix(rep(NA, n*m), nrow=n, ncol=m))
names(responses) <- c(template$name, 'obj')

normVals = responses

### NOTE: this takes around 2 minutes to complete on a good PC. It may needs up to 5 minutes on older PCs.
s = Sys.time()
for ( i in 1:n ) {
  p <- createCandidateParam(template)
  fitness <- getResponse(p, simulation, objectiveFunc)
  responses[i,] = c(p$value, fitness)
  normVals[i,] = c(((p$value-p$min)/(p$max-p$min)), fitness)
}
e = Sys.time()
print (e-s)

# model = clrGetField(simulation, 'tsModel')
# clrReflect(model)
# blah = clrCall(model, "get_KSurf")
# blah

head(responses)

### Let's visualise some results. There is no use of rClr below as such, but this section an important part of the use case.

library(ggplot2)
dev.new()
d <- ggplot(responses, aes_string(x=names(responses)[2], y='obj'))
d + geom_point()

library(reshape)
normVals$.row = rownames(normVals)

resp2 = (melt(normVals[,c(template$name, '.row')], id='.row'))
resp2$obj = normVals$obj # NOTE: this line may not be recommended practice...
pcp <- ggplot(resp2, aes(variable, value, group=.row, color=obj))
pcp + geom_line()

resp2 = resp2[resp2$obj > 0.57,]
pcp <- ggplot(resp2, aes(variable, value, group=.row, color=obj))
pcp <- pcp + ggtitle("Parameters with NSE above 0.57")
pcp + geom_line()

